ELECTROMAGNETIC  ANALYSIS  OF  PERIODIC 
STRUCTURE  REFLECTION  AND  TRANSMISSION 
BY  THE  HYBRID  FINITE  ELEMENT  METHOD 


Daniel  T.  McGrath,  Major,  USAF 


September  1995 


Final  Report 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  IS  UNLIMITED. 


PL-TR-95-1103 


This  final  report  was  prepared  by  the  Phillips  Laboratory,  Kirtland  Air  Force  Base, 

New  Mexico,  under  Job  Order  31 52AK06.  The  Laboratory  Project  Officer-in-Charge  was 
Maj.  Daniel  T.  McGrath  (¥/SR). 

When  Government  drawings,  specifications,  or  other  data  are  used  for  any  purpose  other  than 
in  connection  with  a  definitely  Government-related  procurement,  the  United  States 
Government  incurs  no  responsibility  or  any  obligation  whatsoever.  The  fact  that  the 
Government  may  have  formulated  or  in  any  way  supplied  the  said  drawings,  specifications,  or 
other  data,  is  not  to  be  regarded  by  implication,  or  otherwise  in  any  manner  construed,  as 
licensing  the  holder,  or  any  other  person  or  corporation;  or  as  conveying  any  rights  or 
permissiori  to  manufacture,  use,  or  sell  any  patented  invention  that  may  in  any  way  be  related 
thereto. 


This  report  has  been  authored  by  an  employee  of  the  United  States  Government.  Accordingly, 
the  United  States  Government  retains  a  nonexclusive  royalty-free  license  to  publish  or 
reproduce  the  material  contained  herein,  or  allow  others  to  do  so,  for  the  United  States 
Government  purposes. 

This  report  has  been  reviewed  by  the  Public  Affairs  Office  and  is  releasable  to  the  National 
Technical  Information  Service  (NTIS).  At  NTIS,  it  will  be  available  to  the  general  public, 
including  foreign  nationals. 

If  your  address  has  changed,  if  you  wish  to  be  removed  from  the  mailing  list,  or  if  your 
organization  no  longer  employs  the  addressee,  please  notify  PLAVSR,  3550  Aberdeen 
Ave  SE,  Kirtland  AFB,  NM  87117-5776,  to  help  maintain  a  current  mailing  list. 

This  report  has  been  reviewed  and  is  approved  for  publication. 


DANIEL  T.  McGRATH,  Maj.,  USAF 
Deputy  Chief,  Electromagnetic  Sources  Division 
Project  Officer 


FORREST  J.'  AGEE,  GM-15 
Chief,  Electromagnetic  Sources 
Division 


FOR  THE  command: 


WILLIAM  G.  HECKATHORN,  Colonel,  USAF 
Director,  Advanced  Weapons 
and  Survivability  Directorate 


DO  NOT  RETURN  COPIES  OF  THIS  REPORT  LTsFLESS  CONTRACTUAL  OBLIGATIONS  OR 
NOTICE  ON  A  SPECMC  DOCUMENT  REQUIRES  THAT  IT  BE  RETURNED. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704^0188 


Public  reoomng  burden  for  this  coMeaion  of  information  is  estimated  to  average  1  Hour  per  response,  including  the  time  for  reviewing  instruaions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  colleaion  of  information  Send  comments  regarding  this  burden  estimate  or  any  other  aspea  of  this 
coHeaion  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  HeacJouarters  Services,  Direaorate  tor  information  Operations  and  Repoas,  1215  Jefferson 
Da/ts  Highway.  Suite  1204.  ^'iington,  VR  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduaion  Project  (0704-0188),  Washington,  DC  205C3. 


1.  AGENCY  USE  ONLY  (Leave  blank) 


2.  REPORT  DATE 

September  1995 


3.  REPORT  TYPE  AND  DATES  COVERED 

Final;  Aug  93  -  Sep  95 


4.  TITLE  AND  SUBTITLE 

ELECTROMAGNETIC  ANALYSIS  OF  PERIODIC 
STRUCTURE  REFLECTION  AND  TRANSMISSION  BY 
THE  HYBRID  FINITE  ELEMENT  METHOD 


6.  AUTHOR(S) 

Daniel  T.  McGrath,  Major,  USAF 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADORESS(ES) 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


Phillips  Laboratory 
3550  Aberdeen  Avenue,  SE 
Kirtland  AFB,  NM  87117-5776 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


PL-TR--95-1103 


to.  SPONSORING  /  MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Approved  for  Public  Release; 
Distribution  Unlimited 


13.  ABSTRACT  (Maximum  200  words) 

A  numerical  solution  for  the  time-harmonic  electromagnetic 
fields  in  an  infinite  periodic  structure  has  been  developed,  imple¬ 
mented,  and  validated.  It  uses  a  variation  of  the  hybrid  finite 
element  method  that  includes  periodic  radiation  conditions  at  the 
exterior  surfaces,  and  periodicity  conditions  at  the  unit  cell  side 
walls  inside  the  structure.  The  electric  fields  inside  the  unit 
cell  are  represented  by  linear  vector  finite  elements.  The  success¬ 
ful  test  cases  include  an  inductive  screen,  a  parallel  wire  grid,  a 
conducting- sphere  artificial  dielectric,  a  bandgap  structure  made  up 
of  parallel,  layered-dielectric  rods,  and  crossed  conducting  bars 
representative  of  reinforced  concrete.  The  report  describes  the 
formulation  of  the  method,  and  its  implementation  in  a  general-pur¬ 
pose  computer  code,  and  the  results  of  the  test  cases.  Possible 
uses  of  the  code  for  studies  in  high  power  microwaves  are  addressed. 


14.  SUBJECT  TERMS  15.  PAGES 

Electromagnetic  F.elds,  Scattering,  Filters 

Antenna  Arrays,  Windows,  Filters,  Polarizers  - - - — 

Finite  Element  Analysis,  Computer  Programs 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACT 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 


Unclassified 


Unclassified 


Unclassified 


Standard  Form  298  (Rev  2-89) 

P'escrioed  by  ANSi  Std  Z39-18 
298-102 


NSN  7540-01-280-5500 


TABLE  OF  CONTENTS 


1.  INTRODUCTION  . 1 

2.  BACKGROUND  . . 2 

2.1.  Problem  Description . 2 

2.2.  Solution  Method  Overview  . . 2 

3.  THEORY . .6 

3.1.  The  Weak  Form  of  the  Wave  Equation  . 6 

3.2.  Finite  Element  Representation  of  the  Electric  Field  . . 7 

3.3.  Interior  Discretization . 8 

3.4.  Mode  Functions  and  the  Periodic  Integral  Equation  .  . . 11 

3.5.  Periodic  Radiation  Boundary  Terms  . . 14 

3.6.  Side  Wall  Periodicity  Condition  . . 15 

3.6.1.  Boundary  Wrapping  Concept  . . . . 15 

3.6.2.  Matrix  Transformation  . . . . . . 16 

3.6.3.  Algorithmic  Implementation  . . 17 

3.6.4.  "Overlap  Elements"  on  the  Radiation  Boundaries  ..................  18 

3.7.  Incident  Field  Terms . 18 

3.8.  Transmissivity  and  Reflectivity  . 19 

4.  NUMERICAL  IMPLEMENTATION . 21 

4.1.  Geometry  Generation  . 21 

4.1.1.  CAD  Geometry  Definition . 21 

4.1.2.  Unit  Cell  Side  Wall  Geometry  Definition . 21 

4.1.3.  Conductor  Boundary  Conditions  . 23 

4.1.4.  Node-Based  to  Edge-Based  Geometry  Conversion . 24 

4.2.  Matrix  Data  Structure  and  Computations . 7T  .  .  .  .  25 

4.3.  Matrix  Solution . . . 26 

4.3.1.  Conjugate  Gradient  Method  . 26 

4.3.2.  Biconjugate  Gradient  Method  . . . . .  .  27 

5.  VALIDATION  . 31 

5.1.  Inductive  Screen  Filter  . 31 

5.2.  Wire  Grid  Array . 34 


111 


.  .  .  36 


5.3.  Artificial  Dielectric  .................................... 

5.4.  Dielectric  Bandgap  Structure  .................................  38 

5.5.  Reinforced  Concrete  ......................................  41 

5.6.  Crossed  Dipole  Frequency  Selective  Surface  .......................  43 

5.7.  Storage,  Execution  Time  and  Accuracy  Issues  ......................  45 


6. 


CONCLUSIONS  AND  RECOMMENDATIONS 


48 


REFERENCES 


50 


LIST  OF  FIGURES 

1.  Representative  Problem  Geometry:  (a)  Side  view  showing  upper  and  lower  radiation 

boundaries;  (b)  top  view  of  array  lattice  showing  unit  cell  boundaries  . . 3 

2.  Example  Tetrahedron  Mesh:  Crossed  Dipole  FSS  Element  in  Skewed  lattice  ......  4 

3.  Side  View  of  Example  Mesh  -  Surface  Grid  on  Unit  Cell  Walls  Only  ..........  22 

4.  Direction  and  Sequence  for  Defining  Curves  and  Mesh  Areas  to  Achieve  Identical  Meshes 
on  Opposing  Unit  Cell  Side  Walls:  (a)  Curve  Directions;  (b)  Order  of  Selection  for  Curves 
to  Make  up  Mesh  Areas  ......................................  23 

5.  Convergence  of  Residual  Norm,  Reflectivity  and  Transmissivity  using  Conjugate  Gradient 


14.  Layered  Dielectric  Rod  used  to  Simulate  Photonic  Band  Structure . . . 39 

15.  Tetrahedron  Mesh  for  Layered  Dielectric  Rods:  (a)  Exploded  View  of  Single 

Rod;  (b)  Unit  Cell  Mesh  for  Four-layer  Structure . . . 40 


16.  Comparison  of  Hybrid  Finite  Element  Method  Calculations  and  Measurements 

for  10-layer  Dielectric  Rod  Array . 41 

17.  Cutaway  of  One  Half  of  Tetrahedron  Mesh  for  Reinforced  Concrete  Unit  Cell  ....  43 

18.  Normal-incidence  Transmission  through  4''x4"  Reinforcing  Bar  Lattice  . 44 

19.  Oblique-incidence  Transmission  through  4"x4”  Reinforcing  Bar  Lattice: 

(0o,0o)=(3O°,45°),  Transverse  Electric  Polarization  . . . 44 

20.  Crossed  Dipole  FSS  Element  and  Lattice  Geometry  . 46 

21.  Comparison  of  HFEM  and  MoM  Calculations  of  TM  Reflection  from  Crossed  Dipole 

Frequency  Selective  Surface  (6q=30° ,(f>Q=0°)  . . 46 

22.  Tetrahedron  Mesh  for  Crossed  Dipole  Frequency  Selective  Surface  Test  Case: 


Two  Air  Layers  Surround  the  Conducting  Elements  . . 47 

23.  Number  of  Iterations  vs.  Frequency  for  Bandgap  Structure . 49 


LIST  OF  TABLES 

Table  1.  Mesh  Sizes  and  Storage  Requirements  for  Test  Cases  . 48 

Table  11.  Iterations  and  Execution  Times  for  Test  Cases . 48 


\ 


1 1  I 


V 


1.  INTRODUCTION 


Analysis  and  design  of  periodic  structures  is  important  to  many  areas  of  microwave  and 
antenna  engineering.  For  example,  in  designing  equipment  enclosures  to  suppress  electromag¬ 
netic  interference  (EMI),  one  may  need  to  know  the  amount  of  radio  frequency  (RF)  penetration 
through  periodic  screens  such  as  perforated  metal,  wire  grids,  or  reinforced  concrete.  In  anten¬ 
na  design,  periodic  structures  may  take  the  form  of:  frequency  selective  surfaces  (FSSs)  [1],[2], 
proposed  for  use  as  subreflectors  and  interference-suppressing  radomes;  polarizers  [3];  artificial 
dielectrics  [4]  for  lightweight  lenses;  and  angular  filters  for  sidelobe  suppression  [5].  Periodic 
structures  also  appear  in  the  field  of  optics,  as  diffraction  gratings  and  bandgap  filters.  But  in 
optics,  as  in  radio  frequency  applications,  the  behavior  of  these  devices  is  essentially  electromag¬ 
netic.  In  most  cases,  acceptable  performance  prediction  is  obtainable  only  with  analytical  meth¬ 
ods  that  account  for  mutual  interactions  between  elements  of  the  structure.  The  available  tools 
have,  to  date,  been  based  on  mode  matching  and  moment  methods.  Their  most  significant 
limitations  were:  (1)  inability  to  model  "inhomogeneous"  dielectrics  (anything  other  than  contin¬ 
uous  planar  slabs);  and  (2)  restrictions  on  the  shape  and/or  orientation  of  conducting  elements, 
requiring  a  separate  computer  code  for  each  of  many  classes  of  geometries.  This  report  de¬ 
scribes  a  new  analysis  technique  that  employs  the  hybrid  finite  element  method  (HFEM)  in  order 
to  permit  full  wave  analysis  of  very  general  periodic  structures  containing  dielectric  inhomogene¬ 
ities  and  arbitrarily-oriented  conductors.  It  is  an  extension  of  previous  work  in  modeling  phased 
array  antennas  [6]. 
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2o  BACKGROUND 


The  generic  problem  to  be  addressed  is  electromagnetic  transmission  through  large, 
planar  periodic  structures.  They  are  assumed  to  be  large  enough  that  edge  effects  are  negligible, 
allowing  use  of  an  infinite  structure  model  This  restricts  the  analysis  to  a  single  unit  cell  that 
is  typically  a  fraction  of  a  wavelength  on  each  side  at  the  frequencies  at  which  the  device  is 
designed  to  operate.  Figure  1  illustrates  the  general  problem  geometry  and  defines  the  notation 
conventions  used  throughout  this  report. 

The  structure  may  include  conducting  obstacles  with  arbitrary  shape,  and  any  number  of 
dielectric  regions  with  distinct  permittivity,  e,  and  permeability,  p..  The  dielectrics  are  isotropic 
(their  properties  are  not  a  function  of  direction)  and  linear  (their  properties  do  not  depend  on  the 
magnitude  of  the  fields).  The  structure  lattice  may  be  skewed,  as  illustrated  in  Figure  lb. 


2o2o  Solutioia  Method  Overview 

The  solution  technique  is  a  hybrid  of  the  finite  element  method  (FEM)  with  an  integral 
equation  representation  for  fields  above  and  below  the  structure.  The  unit  cell  is  truncated  at 
planes  of  constant  z  above  and  below  the  structure  so  that  all  of  the  materials  (dielectrics  and 
conductors)  are  between  the  surfaces  denoted  T^.  and  ,  which  will  be  referred  to  as  the 
lower  and  upper  radiation  boundaries,  respectively.  The  volume  thus  formed  by  the  truncated 
unit  cell  will  be  denoted  by  0. 

The  free  space  and  dielectrics  inside  0  will  be  subdivided  into  small  volume  elements 
(tetrahedra),  referred  to  as  "cells."  The  electric  field  will  be  represented  in  terms  of  simple 
expansion  functions  or  "finite  elements"  defined  relative  to  cells.  The  particular  form  used  in 
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Figure  1.  Representative  Problem  Geometry:  (a)  Side  view  showing  upper  and 
lower  radiation  boundaries;  (b)  top  view  of  array  lattice  showing  unit  cell  boundaries 
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perfectly  conducting  surfaces.  Assuming  that  the  edges  are  ordered  so  that  the  first  are  in 
and  the  last  are  in  F^^  ,  then  the  matrices  S^‘  and  S^'*'  due  to  the  periodic  radiation  conditions, 
are  zero  except  for  dense  upper  left  and  lower  right  submatrices,  respectively.  Finally,  multipli¬ 
cation  on  the  left  by  the  matrix  R  and  on  the  right  by  its  Hermitian  (conjugate  transpose)  imple¬ 
ments  the  side  wall  periodicity  condition.  The  right  side  vector  is  due  to  a  unit-amplitude 
plane  wave  (the  dominant  Floquet  mode).  Its  entries  are  zero  except  for  those  corresponding 
to  edges  in  F^. .  Solution  of  the  system  of  equations  gives  the  vector  E,  from  which  reflectivity 
and  transmissivity  are  computed. 

The  following  chapter  will  give  the  details  of  how  each  of  the  terms  in  the  above  matrix 
equation  is  derived  from  electromagnetic  theory  and  from  the  methods  of  finite  elements. 
Chapter  4  will  discuss  how  it  was  implemented  in  a  general-purpose  computer  code.  Chapter 
5  will  show  the  results  for  several  validation  and  demonstration  cases. 
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So  THEORY 

3olo  The  Weak  Form  of  the  Wave  Equation 

Finite  element  solutions  are  usually  based  on  discretization  of  a  "variational  principle" 
or  of  a  "wealc  form."  Application  of  the  Rayleigh-Ritz  method  to  the  variational  principle 
results  in  the  same  matrix  equation  as  does  application  of  GalerMn’s  method  to  the  weak  form 
[6:102-103]. 

The  wave  equation  for  the  electric  field  in  linear,  isotropic,  inhomogeneous  media  is 


VxlVxE-kle.E  =  0 

IJ'r 


(2) 


where  kg  is  the  free  space  wavenumber,  and  the  relative  permeability,  and  permittivity,  e^, 
may  be  functions  of  position.  The  inner  product  of  (2)  with  a  weighting  function  W  is 


Vx±VxE-kle^E 


dv  =  0 


(3) 


Using  a  Green’s  identity: 


J_Vx « Vx 
k-r 


E- 


J_W*  xVxEo  fids  =  0 
k-r 


(4) 


where  F  is  the  union  of  all  boundaries  enclosing  0  as  well  as  conducting  surfaces  inside  0. 
Equation  (4)  is  called  a  "weak  form"  because  the  Green’s  identity  shifted  one  of  the  derivatives 
from  E  to  W,  weakening  the  differentiability  requirement  on  E.  The  final  form  is  obtained  by 
using  the  Maxwell  curl  equation  VxE  =  and  the  constitutive  relation  wpt  =  kgr}  : 


6 


I 


—Vxw*  -VxE-kle^W* 

l^r 


E 


dv-jkQVo\  W*  ‘  (nxH)ds  =0 
r 


(5) 


The  term  on  the  right,  a  surface  integral  over  all  boundaries  enclosing  0,  provides  the  mecha¬ 
nism  for  enforcing  boundary  conditions.  On  the  open  parts  of  F,  suitable  expressions  for  the 
tangential  magnetic  field  may  be  substituted  to  impose  radiation  boundary  conditions,  as  dis¬ 
cussed  in  Section  3.4. 

For  perfect  conductors,  the  second  integral  above  becomes  zero.  The  tangential  electric 
field  is  zero,  so  admissible  functions  for  E,  and  hence  also  for  W,  are  those  whose  tangential 
component  vanishes  at  a  perfectly  conducting  boundary. 


3.2.  Finite  Element  Representation  of  the  Electric  Field 

The  electric  field  will  be  expanded  in  "finite  elements,"  which  are  functions  defined 
locally  within  each  tetrahedral  cell.  The  "linear  edge  elements"  are  defined  relative  to  two  mesh 
nodes  i  and  j  that  bound  a  mesh  edge  [7]-[9]: 

(6) 

where  L^j  is  the  length  of  the  edge.  The  scalar  functions  and  fj  have  the  value  1  at  node  i  or 
j,  respectively,  and  decrease  linearly  to  0  ar  the  other  three  nodes  in  the  cell.  The  vector  func¬ 
tion  is  normal  to  the  faces  opposite  nodes  i  and  j.  Its  component  along  edge  ij  is  a  unit 
vector.  Now  the  electric  field  is  approximated  by 

£  »  E  =  (7) 

5=1 

where  is  a  scalar,  complex  coefficient  representing  the  field  magnitude  and  phase  along  the 
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edge  s.  This  expansion  makes  it  a  trivial  matter  to  enforce  boundary  conditions  along  perfect 
conductors,  even  at  sharp  edges  and  points.  It  also  ensures  that  spurious  modes  will  not  be 
generated  [10].  Furthermore,  it  correctly  enforces  continuity  of  the  normal  electric  field  across 
dielectric  interfaces  [1 1].  All  three  of  these  are  reasons  it  is  preferred  over  the  scalar  expansion 
E  «  (defined  relative  to  nodes). 


The  elements  of  the  sparse  matrix  in  (1)  are  due  to  the  volume  integral  term  in  (5). 
Substituting  the  expansion  (7)  for  E  and  using  each  in  turn  as  a  weighting  function  (Galerkin’s 
method): 


(8) 


where  0^^  denotes  the  union  of  all  cells  that  share  edges  s  and  t. 

The  integrations  in  (8)  are  carried  out  with  the  help  of  a  coordinate  transform  local  to 

each  tetrahedron  (called  "simplex"  or  "barycentric"  coordinates  )  [12:266-274].  Let  U  denote 

the  matrix  formed  using  the  four  vertex  coordinates: 

r  1 


u  = 


1 

1 

i 


Xi 

X2 

X3 

X4 


3^1 

Jl  ^2 

}'3  23 

)'4  24 


(9) 


and  let  T  be  the  matrix  of  the  16  cofactors  of  U.  There  are  four  coordinate  directions  in 
the  tetrahedron  given  by 
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(10) 


where  V  is  the  cell  volume.  Since  r,-  is  equal  to  unity  at  node  i  and  goes  linearly  to  zero  at  all 
other  nodes,  it  is  the  correct  form  for  the  scalar  finite  element  function,  and  .  The 

gradient  of^-  is 

The  curl  of  the  finite  element,  V  x  is  also  needed  in  evaluating  the  integrals  in  (8).  For  that, 
the  following  vector  identity  is  used: 

(12) 

Note  that  the  first  term  on  the  right  in  (12)  is  zero  since  it  is  a  second  derivative  of  a  linear 
function.  Consequently, 

vx?:,  =  vx(/;.v^-/jvy;.) 

-2V/jXV^. 

To  evaluate  the  contribution  of  a  cell  k  to  (8),  let  i  and  j  denote  the  local  indices  of  the 
nodes  bounding  edge  s;  and  let  m  and  n  denote  those  for  edge  t.  All  four  of  these  local  node 
indices  are  between  1  and  4.  The  first  integral  term  in  (8)  will  be 

s'm  -  X  V/^  •  V/„  X  V/„  (14) 

•  [  (^i35}4  -  T,4Tp)iJ„,T„^  -  ) 

+  {TiJj2-Ti2Tj^){T^Jn2-Tm2TnA)  (15) 

(^m2^n3  ~  '^m3^n2)  ] 
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And  the  second  volume  integral  term  in  (8)  will  be 


Qjt 

-fifn^  “  ^fm  ^  fjfn^fi  ”  ^fm 


(16) 


Since  the  gradient  terms  in  (15)  are  constants,  they  may  be  taken  outside  the  integral,  leaving 
four  terms  of  the  form 


hj  =  [  fifjdxdydz 

k 

which  transform  to  simplex  coordinates  as  [12] 


V 


j  tj  dtyfj 
0 


(17) 


(18) 


where  1/6!^^  is  the  Jacobian  of  the  transformation.  This  integral  evaluates  to 

where  6,---  is  the  Kronecker  delta.  Finally,  (16)  reduces  to 

V 

S'ft  =  - E  [ ( 1  - «/») Tj,T„,  T„T„  ,20) 

Each  matrix  entry  is  the  sum  of  contributions  from  (15)  and  (20)  from  all  cells  that  share 
edges  s  and  t.  Thus,  those  two  equations  give  the  interior  matrix  entries  as  closed-form  algebra¬ 
ic  expressions  that  depend  only  on  the  cell  geometry,  its  constitutive  parameters,  and  the  wave- 
number.  The  following  two  sections  show  the  derivation  for  the  entries  in  the  matrices  pertain- 
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ing  to  the  radiation  boundaries. 


3.4.  Mode  Functions  and  the  Periodic  Integral  Equation 

The  top  and  bottom  of  the  unit  cell  are  free  space  boundaries  where  a  radiation  condition 
must  be  applied.  The  method  for  doing  so  is  to  substitute  an  integral  equation  for  nxH  into  the 
boundary  term  of  the  weak  form.  Due  to  the  periodic  nature  of  the  problem,  the  integral 
equation  may  be  written  as  a  summation  over  "Floquet  modes"  [13]. 

Let  the  electric  field  above  and  below  the  structure  be  expanded  in  a  set  of  mode  func¬ 
tions  where  the  subscript  p  is  1  for  transverse  electric  (i.e.  Elz)  or  2  for  transverse 
magnetic  (i.e.  Hit): 


T  J  *  >^ymny) 

^pmn  ^ 


(21) 


(22) 


(23) 


^xmn  =  hsmOoCOS  <l>0 


Kmn‘k^ sin  sin 2xmcoty 


where  (0o»^o)  the  angle  toward  which  the  incident  wave  is  propagating,  {m,n)  is  the  mode 
index  and  y  is  the  lattice  skew  angle. 
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At  any  point  above  the  structure,  the  transverse  electric  field  may  be  expanded  in  modes 


traveling  in  the  +z  direction: 


2  GO  CO 

£'(z>o)  =  E  E  E  ^pmn  Spmn  ^ 

/?  =  !  m  =  -oo  n  =  -oo 


‘J^mn  ^ 


(26) 


-k^  - 

^0  ^xmn  ymn 


1/2 


(27) 


C^mn  is  an 


coefficient.  The  transverse  magnetic  field  in  each  mode  is 


pmn  pmn 


(28) 


where  is  a  modal  admittance: 


Y 


Imn 


^mn 

kgrio 


(29) 


2mn 


kn 


%n  ^0 


(30) 


The  total  magnetic  field  above  the  structure  is 


Z  XH'(z>0)  =  -E  E  E  c;„pypmngpn,y''‘’"'  (31) 

pmn 

In  order  to  write  the  right  side  in  terms  of  E  \  the  orthogonality  property 

J  Spmn  °  ^qij  ~  ^ pq^mi^ nj  (32) 

unit  cell 

is  used  to  solve  (26)  for  the  coefficients: 
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^pmn  I  ^  ^Spmn^^y  (33) 

unit  cel! 

with  the  result 

Z-XH'(Z>0)  =  -E  E  E  yp^nSpn..  I  E‘-g;„„<bdy  (34) 

^  ^  ”  unit  cell 

On  the  lower  unit  cell  boundary,  the  transverse  electric  field  includes  an  incident  field 
term  in  one  of  the  two  lowest  order  modes  depending  on  whether  it  is  transverse  electric  {q=  1) 
or  transverse  magnetic  (q=2): 

E'u<o) - -i:  i:  E  os) 

p  =  l  m  =  -oo  /i  =  -oo 

The  signs  of  the  exponential  terms  indicate  the  direction  of  propagation.  Again,  using  the 
modes’  orthogonality  property, 

<^pmn=  I  -Spmn^dy-Sp^S^o^nO  (36) 

unit  cell 

giving  the  magnetic  field  below  the  structure  as 

zxH'(z<0)  =  -E  E  E  JV™.  «/>-»«  f  E'-g;„p<lxdy  e-"-"" 

P  "  "  uJicai  (37) 

*  2  l',ooi,oo 

The  above  expressions  (34)  and  (37)  are  the  integral  equations  that  are  substituted  into  the  weak 
form  (5)  to  enforce  the  radiation  conditions.  The  following  section  derives  expressions  for  the 
resulting  matrix  terms. 
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The  matrix  entries  pertaining  to  the  upper  radiation  boundary,  are  found  by  substi¬ 
tuting  ^  for  in  (34),  then  substituting  (34)  for  ny-H  and  for  W*  in  Eq.  (5): 


^st  ~  SEE  ^pmnSpmn^^y  [  °  Spmn}^  ^ 

J  p  m  n  ^ 


(38) 


'  z-^ 


&ch  integral  term  contains  a  term  in  exp{i;j(^^^  +  ky^^)}  ,  which,  besides  the  finite  ele¬ 
ments,  is  the  only  term  with  x  or  y  dependence.  Let  ^({k^,ky)  denote  the  two-dimensional 
Fourier  transform  of  that  part  of  lying  on  the  radiation  boundary: 


Spmn^dy  =  ^pmn  “  | 


t  Spmn 


J  (^xmn^  ^  ^ymny  1 


dxdy 


z* 


^pmn  °  ^ t^^xmn’^ymn'^ 


pm 

=  h 


pmn  ^tmn 


(39) 


Now  the  matrix  entries  may  be  written  in  the  more  compact  form 


III  ■  I  ]  •  i 

pmn 


* 

smn 


(40) 


Finally,  the  term  in  brackets  may  be  simplified  by  summing  TE  and  TM  modes: 


E  E  ^tmn  “  ^mn  “  ^ smn 

m  n 


(41) 


,  2  _  ,  2  _  ,2 

^0  ^xmn  ^ymn 


-1/2 


mn 


d-xdykQ  7\q 


(^0  ^ymn^  ^xmn^ymn 
^xmn^ymn  (^0~^xmn^ 


(42) 


The  expression  for  the  matrix  pertaining  to  the  lower  radiation  boundary  is  identical  to  (41). 


This  follows  from  the  fact  that  the  summations  in  (35)  and  (37)  are  the  same  except  for  the  sign 
of  the  exponent  inside  the  integral.  In  evaluating  those  integrals,  the  value  of  z  in  the  exponent 
is  taken  as  zero  because  it  represents  only  a  constant  phase  factor. 

This  concludes  the  derivations  for  the  entries  in  various  parts  of  the  matrix  S.  The 
following  section  discusses  the  matrix  transformation  that  imposes  a  periodicity  condition  on  the 
open  unit  cell  side  walls. 

3.6.  Side  WaU  Periodicity  Condition 
3.6.1.  Boundary  Wrapping  Concent 

To  understand  how  the  transformation  matrix  R  in  (1)  is  created,  it  is  important  to  recog¬ 
nize  that  the  unit  cell  mesh  is  actually  part  of  an  "infinite  mesh."  A  portion  of  the  mesh  is 
replicated  identically  within  each  unit  cell.  An  essential  requirement  is  that  the  mesh  continues 
unbroken  across  unit  cell  boundaries,  which  implies  that  the  mesh  on  a  unit  cell  side  wall  is  an 
exact  replica  of  that  on  the  opposing  side  wall.  Section  4.1.2  will  discuss  how  to  accomplish 
this  "mesh  continuity"  using  a  typical  commercial  mesh  generation  program. 

The  system  of  equations  formed  by  (15),  (20)  and  (41)  contains  one  unknown  per  mesh 
edge.  But  because  of  the  Floquet  condition,  the  fields  on  opposite  side  walls  are  identical  except 
for  a  constant  and  known  phase  shift  that  depends  only  on  the  unit  cell  size  and  the  scan  angle. 
Consequently,  the  linear  system  is  overdetermined  because  there  is  a  linear  dependence  between 
unknowns  associated  with  image  edges.  For  example,  if  e^  is  the  electric  field  along  an  edge 
in  then  the  field  along  the  corresponding  edge  in  F^.  is 
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(43) 


„  c  in  ( 


e^/  =  e^e"  "  *  ,  sE  ,  s  E 


/3^  =  ^0  sini^o  cos(^o 


(44) 


Similarly,  the  fields  for  corresponding  edges  in  Fy^  and  Fy.  are  related  by 

^  ,sEVy^,  s^EFy. 


(45) 


^y  =  kQsindQ  sm<f)Q  (46) 

Because  of  this  dependence,  the  linear  system  may  not  be  solved  until  it  is  reduced  by  combining 
dependent  rows  and  columns. 

One  way  of  viewing  the  reduction  procedure  is  the  following:  The  infinite  mesh  will  be 
recreated  from  the  unit  cell  mesh  by  "wrapping"  opposite  sides  onto  each  other  with  an  approp¬ 
riate  phase  shift.  This  will  merge  the  image  edges,  thus  eliminating  dependent  unknowns  from 
the  linear  system.  There  are  two  implementations:  the  matrix  transformation  indicated  in  Eq. 
(1),  and  algorithmic  implementation. 

3.6.2.  Matrix  Transformation 

The  matrix  transformation  is  executed  by  multiplying  on  the  left  by  a  reduction  matrix 
R  and  on  the  right  by  its  Hermitian.  This  relies  on  a  special  ordering  of  the  unknowns,  with 
those  in  the  side  walls  last,  so  that  they  make  up  the  lower  right  portion  of  the  matrix  S. 

R  is  MxN,  where  N-M  is  the  number  of  edges  on  F^^  and  Fy^  combined.  Its  first  M 
columns  are  the  MxM  identity  matrix.  In  the  remaining  columns,  there  is  an  entry  R^-^  for  each 
pair  of  image  edges; 
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b.  SET  Sg.(.—  Sjij/  +  Sjj 

c.  SET  S^t  =  0 

2.  ELSE  IF  t  IS  NOT  ON  THE  +x  BOUNDARY,  THEN: 

a.  SET  S,,j  =  Sjt  exp{j/?J 

b.  SET  S^,  =  S,j  exp{-j^^} 

c.  SET  S^t  =  S^  =  0 
n.  REPEAT  I  FOR  +y  AND  -y  BOUNDARIES 

m.  COMPRESS  THE  MATRIX  (ELIMINATE  ZERO  ROWS  &  COLUMNS) 
IV.  COMPRESS  THE  INCIDENT  FIELD  VECTOR 
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3.6.4.  "Overlap  Elements"  on  the  Radiation  Boundaries 


The  two  periodicity  conditions,  one  for  side  walls  and  the  other  for  radiation  boundaries, 
must  work  in  concert,  neither  duplicating  nor  opposing  each  other.  Gedney  first  establish  a 
straightforward  method  for  doing  so  in  two  dimensions,  which  he  called  "overlap  elements" 
[14].  In  computing  the  terms  of  (1)  pertaining  to  the  radiation  boundaries,  those  edges  in  rx^ 
and  Fy^.  are  ignored.  There  are  no  terms  in  or  S^'  for  those  edges.  However,  edges  in 
Fj^.  and  Fy.  are  treated  as  though  they  were  part  of  triangular  elements  extending  into  adjacent 
unit  cells.  That  is,  the  finite  elements  connected  to  those  edges  "overlap"  the  unit  cell  bound¬ 
aries.  This  method  actually  results  in  a  simpler  algorithm  than  Eq.  (1)  suggests,  since  the  matrix 
transform  need  only  be  applied  to  S^. 

3o7o  Imcndent  FseM  Terms 

The  right  side  vector  in  the  matrix  equation  arises  from  the  second  term  on  the  right  side 
of  (37),  due  to  a  unit  amplitude  incident  electric  field  arriving  from  z<0.  Its  form  is  a  plane 
wave,  which  is  the  0,0-order  Floquet  mode.  Along  edge  t: 

~  I  ^t°  ^ qQoS  qQo)  dX  dy 

r 

‘■z- 

where,  again,  the  constant  phase  factor  exp{-y  z]  is  ignored.  Taking  all  of  the  terms  outside 
the  integral  that  do  not  depend  on  x  and  y: 

inc  ,7  I  f  T  j^^xOO^*^yOQy'>  ,  , 

Cf  =  "-2;^or;oF^oo  Vo  °  J  (49) 

Tz- 

Note  that  the  integral  is  the  two-dimensional  Fourier  transform  of  the  finite  element’s  tangential 
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component  on  F^.,  written  as  ^(Qq.  Note  also  that: 


sin^^o 


(50) 


cos  00 

Vo 

1 

Vo  cos  00 


q  =  l 
q  =  2 


and  now  (49)  may  be  rewritten  in  the  form 


(51) 


cos0o(isin<^o  ■)’cos<^)o)  •  1,00  ^  =  1 

sec 00 (i cos 4>q  +y sin  <^o )  ’  i^roo  9  =  2 


(52) 


The  solution  of  the  resulting  system  of  equations  gives  the  vector  of  unknown  coefficients  E  in 
(1).  The  following  section  shows  how  these  edge-field  coefficients  are  used  to  find  the  transmis¬ 
sion  and  reflection  coefficients  for  the  structure. 


3.8.  Transmissivity  and  Reflectivity 

The  most  important  quantities  to  be  found  from  a  periodic  structure  analysis  are  the 
surface  transmissivity  and  surface  reflectivity.  These  complex  quantities  provide  knowledge  of 
how  the  incident  power  divides  itself  between  transmitted  and  reflected  components,  as  well  as 
the  phase  delay  of  the  transmitted  component. 

The  transmission  coefficients  for  each  mode  are  given  by  (33),  since  C^mn  represents 
the  excitation  of  each  mode  propagating  away  from  the  array  in  the  +z  direction.  However, 
when  the  transmission  is  into  a  different  mode  than  the  incident  wave,  such  as  in  the  case  of 
grating  lobes  and  Bragg  lobes,  the  transmission  coefficients  need  to  take  into  account  the  differ- 
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ence  in  modal  admittances.  Letting  denote  the  transmission  coefficient: 


pmn 


A 


Y. 


P”'”  p”  .h 
Y  ^ smn  pmn 

'«00 


(53) 


The  greatest  interest  is  usually  in  the  lowest  order  transmission  coefficient  in  the  m=n=0  mode, 
especially  when  d^,dy<.5\,  in  which  case  all  higher  order  Floquet  modes  are  evanescent.  The 
lowest  order  transmission  coefficient  is  the  transmissivity,  given  by 


5j00 


1 


Y. 


poo 


^qOO 


(54) 


X  sin  <l)Q  -  y  cos  (t>Q  p  =  1  (55) 

X  cos  4>q -i- y  sin  4>q  p  =  2 

The  modal  reflection  coefficients  are  found  from  (36).  Except  for  the  m=n=0  mode,  they  are 

as  given  by  (48)  above,  with  F-  replacing  Fd-.  The  lowest  order  coefficient,  the  reflectivity,  is 


^pOO 


^sOO  ‘  Pr,-^ 


Pq 


Y. 


poo 


(56) 


This  concludes  the  presentation  of  the  theory  and  formulation.  The  next  chapter  address¬ 
es  various  aspects  of  the  implementation  in  a  general-purpose  computer  code. 
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4.  NUMERICAL  IMPLEMENTATION 


4.1.  Geometry  Generation 

4.1.1.  CAD  Geometry  Definition 

The  finite  element  geometries  are  generated  using  a  commercial  "CAD"  software  package 
such  as  I-DEAS™  [15].  The  objective  is  to  create  a  numerical  description  of  a  complex  geome¬ 
try  that  the  electromagnetic  analysis  code  can  interpret  unambiguously.  The  general  steps  the 
user  takes  in  the  geometry  generation  process  are: 

1.  Construct  a  wireframe  geometry  of  points,  lines  and  arcs 

2.  Use  the  curves  (lines  and  arcs)  to  define  surfaces 

3.  Group  curves  into  mesh  areas 

4.  Group  mesh  areas  into  mesh  volumes 

5.  Generate  the  tetrahedral  elements  by  subdividing  mesh  volumes 

6.  Tag  elements  (cells)  with  material  properties 

7.  Tag  nodes  with  boundary  condition  flags 

Steps  1-5  are  fairly  generic,  although  the  details  may  be  peculiar  to  the  particular  software 
package  used.  Steps  6  and  7  will  be  described  in  further  detail  in  sections  4,1.3  and  4.1.4. 

4.1.2.  Unit  Cell  Side  Wall  Geometry  Definition 

A  troublesome  part  of  the  periodic  structure  geometry  generation  is  ensuring  that  the 
mesh  on  opposite  unit  cell  faces  lines  up  correctly.  In  order  for  the  side  wall  periodicity  condi¬ 
tion  to  work,  each  side  wall  edge  must  have  an  identical  "image"  on  the  opposite  side  wall,  with 
the  same  length  and  orientation.  Each  edge  on  the  -hx  boundary  must  have  an  image  on  the  -x 
boundary  that  has  the  same  node  coordinates  except  for  a  translation  of  (<4,0,0).  Each  edge  on 
the  +y  boundary  must  have  an  image  on  the  -y  boundary  that  has  the  same  node  coordinates 
except  for  a  translation  of  (<4cot7,<4,0).  Figure  3  shows  an  example  in  which  all  of  the  geome¬ 
try  is  blanked  out  except  for  the  mesh  on  the  unit  cell  perimeter,  which  is  viewed  from  a  point 


21 


before  the  top  and  bottom  or  any  parts  of  the  interior.  This  ensures  that  those  areas  will  be 
meshed  first,  so  the  size  and  orientation  of  their  triangular  elements  will  not  depend  on  interior 
details.  Furthermore,  the  ordering  of  these  mesh  areas  should  be  such  that  opposite  walls  are 
meshed  sequentially,  for  example,  -x,  +x,  -y,  and  +y.  Finally,  the  curves  making  up  opposite 
side  wall  mesh  areas  must  be  grouped  in  a  consistent  order.  Figure  4b  illustrates  this  by  show¬ 
ing  the  order  that  curves  are  selected  to  make  up  the  +x  and  -x  side  walls.  These  steps  are 
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Figure  4.  Direction  and  Sequence  for  Defining  Curves  and  Mesh  Areas  to  Achieve  Identical 
Meshes  on  Opposing  Unit  Cell  Side  Walls:  (a)  Curve  Directions;  (b)  Order  of  Selection 

for  Curves  to  Make  up  Mesh  Areas 

generally  adequate  to  ensure  that  I-DEAS™  will  produce  identical  grids  on  opposite  unit  cell 
side  walls. 

4.1.3.  Conductor  Boundary  Conditions 

There  are  various  ways  that  any  CAD  package  may  be  used  to  tag  nodes  in  order  to 
signal  the  electromagnetic  code  that  they  represent  conductors.  One  easily-overlooked  difficulty 
pertains  to  the  way  conducting  boundaries  are  represented  by  the  edge-based  finite  elements. 
If  a  mesh  edge  lies  along  a  conductor,  its  tangential  electric  field  must  be  zero.  Therefore,  there 
will  not  be  any  entries  in  the  system  of  equations  or  in  the  matrix  for  those  edges.  However, 
there  can  be  a  field  along  an  edge  that  joins  two  parts  of  the  same  conductor  through  free  space 
or  dielectric.  Consequently,  it  is  not  adequate  to  simply  tag  all  conductor  nodes  the  same.  The 
approach  adopted  in  this  work  is  to  allow  each  node  to  belong  to  any  or  all  of  a  range  of  con- 
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ducting  surfaces,  with  a  different  tag  for  each.  For  example,  the  crossed  dipole  element  in 
Figure  2  consists  of  two  conductor  groups;  one  for  the  horizontal  arms;  and  another  for  the 
vertical  arms.  The  nodes  at  the  junction  will  belong  to  both  groups.  This  allows  the  existence 
of  a  field  component  along  diagonal  edges  joining  vertical  and  horizontal  arms  n^  the  junction. 

4.1.4.  Node-Based  to  Edge-Based  Geometry  Conversion 

Commercial  CAD  programs,  usually  geared  toward  mechanical  engineering,  generate 
geometry  descriptions  in  terms  of  mesh  nodes.  Their  typical  output  consists  of  a  listing  of  node 
coordinates,  then  a  listing  of  the  node  indices  that  define  each  element.  The  edge-based  geome¬ 
try  adds  a  listing  of  edges  by  pairs  of  nodes.  Another  useful  data  structure  is  a  listing  of  ele¬ 
ments  and  the  six  edges  that  comprise  each. 

It  is  convenient,  for  this  periodic  structure  problem,  to  order  the  edges  so  that  those  in 
the  two  radiation  boundaries  are  first  and  last.  As  discussed  in  Section  2.2,  this  causes  the  two 
matrices  associated  with  those  boundaries  to  be  zero  except  for  dense  upper  left  or  lower  right 
submatrices. 

A  crucial  consideration  is  that  the  vector  edge  elements  have  direction.  Throughout  the 
matrix  computations  those  directions  must  be  accounted  for.  Referring  to  Eq.  (6),  the  vector 
function  is  directed  from  the  node  i  to  the  node  j.  A  useful  convention  is  to  define  the  direction 
of  an  edge  as  being  from  the  lower  to  the  higher  global  node  index. 

While  the  matrix  solution  finds  the  field  values  along  mesh  edges,  it  is  often  useful  to 
calculate  the  field  at  nodes.  Node  values  are  compatible  with  the  commercial  CAD  software, 
which  can  be  used  to  generate  "rendered"  plots  of  field  lines  or  field  contours  superimposed  on 
the  geometry.  The  conversion  procedure  amounts  to  reconciling  the  two  approximations 
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(57) 


Ne 

^  *  E  *  E  ^rft 

s-l  t=l 

where  is  the  conventional  node-based  finite  element. 

Consider  the  field  within  a  single  tetrahedron  at  a  coordinate  r.  If  7  is  much  closer  to 
node  i  than  to  any  of  the  other  three  nodes  j,  k,  and  I,  the  value  of  is  nearly  unity,  while  fj, 
and  fi  are  nearly  zero.  I'hen  the  nodal  approximation  to  the  field  is 

E(r)  =  Cl  ft  (58) 

The  field  due  to  the  edge  expansion  is  almost  entirely  due  to  the  three  edges  leading  into  node 
i  since  that  of  the  three  opposing  edges  depends  only  on  fj,  and  ff. 

Efr)  =  E,J  eijf  Vf  .  Ct^f  Vf„  *  i,.,  et,f  Vf  (59) 

Equating  (58)  and  (59)  gives 

‘i  =  EtjetjVfj  *Lt,e,„Vf  *Lt,e„'^f  (60) 

Due  to  the  piecewise  linear  nature  of  the  approximation  for  the  electric  field,  the  right  hand  side 
of  (60)  may  be  slightly  different  within  each  cell  bordering  on  the  node,  so  the  conversion 
procedure  takes  the  average  over  all  cells  adjacent  to  the  node. 

4.2.  Matrix  Data  Structure  and  Computations 

To  take  advantage  of  the  matrix  sparsity,  an  iterative  solution  is  anticipated.  This  is  as 
much  for  reasons  of  storage  as  for  solve  time.  An  enormous  savings  in  computer  memory  can 
be  realized  by  storing  only  the  nonzero  matrix  elements.  However,  a  "direct"  solver,  based  on 
Gaussian  elimination  (usually  LU  decomposition)  creates  new  elements  as  the  solution  proceeds, 
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so  that  the  LU  factors  are  less  sparse  than  the  original  matrix.  Iterative  solvers  such  as  the 
conjugate  gradient  method  (CGM)  [16]  and  biconjugate  gradient  method  (BCGM)  [17],  discussed 
further  in  Section  4.3,  operate  by  making  successive  guesses  at  the  solution  vector,  always  using 


the  original 

execution  ti 
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stored  as  a ' 

each  entry. 

one-dimensional  array,  with  two  integer  arrays  giving  the  row  and  column  index  of 

Of  course,  there  are  sparse  matrix  storage  schemes  that  get  by  with  only  one  integer 

array  [181. 

However,  it  is  most  convenient  to  execute  the  finite  element  computations  one  cell 
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vector  Rq  =  B  -  AXq.  The  initial  search  direction  and  gradient  vectors  Pq  and  Gq  are 

/>„  =  G„  =  (61) 

(The  superscript  H  denotes  Hermitian,  that  is,  conjugate  transpose.)  The  iteration  consists  of 
the  following  steps  using  ||  •  ||  to  denote  the  1^  norm; 


(62) 
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(63) 

(64) 
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A  typical  convergence  test  is  to  examine  the  ratio  e  =  ( ||  i?,- 1|  ^  /  ||  Rq  ||  at  each  iteration, 
A  value  of  10'^  or  10"^  is  usually  adequate  to  ensure  convergence  of  the  transmissivity  and 
reflectivity.  Figure  5  shows  the  convergence  of  the  reflectivity  from  a  single  layeFoT  spherical 
conductors  in  a  square  lattice  (discussed  later  in  Section  5,3),  along  with  the  convergence  of  the 
residual  norm.  Both  of  these  quantities  are  converging  monotonically. 

4,3.2,  Biconjugate  Gradient  Method 

The  generalized  biconjugate  gradient  method  requires  an  additional  two  IxN  complex 
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arrays  for  storage,  which  is  usually  tolerable.  The  search  algorithm  begins  with  the  following 
initial  conditions  [17]: 


M 

o 

(68) 

Rq  -  Pq  -  B 

(69) 

Wq  =  Qq= 

(70) 

The  vectors  W  and  Q  are  in  addition  to  R  and  P,  which  were  required  by  the  CGM.  The 
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following  steps  comprise  the  iteration  procedure: 


(71) 

(72) 

^,>1  =  Ri-OLiAPi 

(73) 

(74) 

Pm  ' 

(75) 

Wm  =  a.i 

(76) 

^•.1  =  (77) 

where  ( • ,  ♦ )  denotes  the  inner  product  of  two  vectors.  Figure  6  shows  the  convergence  of  the 
residual  norm  for  the  same  test  case  as  Figure  5.  Unlike  CGM,  the  residual  norm  does  not 
converge  monotonically,  or  even  gradually.  Nonetheless,  BCGM  required  less  than  .2N  (N  is 
the  number  of  unknowns)  iterations  to  converge  in  this  case,  whereas  CGM  required  1.8N. 
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5.  VALroATION 

The  theory  presented  in  Chapter  3,  and  the  code  design  presented  in  Chapter  4,  were 
validated  by  comparing  the  code  results  for  several  test  cases  to  analytical  formulas,  measure¬ 
ments,  or  calculations  made  with  other,  fundamentally  different  methods.  These  results  serve 
to  show  that  the  new  technique  is  accurate,  flexible  and  dependable.  These  test  cases  were:  (1) 
a  dual-layer  inductive  screen  filter;  (2)  a  wire  grid  array;  (3)  an  artificial  dielectric  made  of 
conducting  spheres;  (4)  a  band-stop  filter  made  of  layered  dielectric  rods;  (5)  a  lattice  of  crossed 
conductor  bars  in  a  dielectric  slab,  emulating  reinforced  concrete;  and  (6)  a  crossed-dipole  FSS. 
The  last  section  in  this  Chapter  compares  the  computational  effort  required  in  terms  of  both 
storage  and  execution  time. 

5.1.  Inductive  Screen  Filter 

The  crossed  conducting  grid  shown  in  Figure  7  is  an  inductive  screen.  (Its  electromag¬ 
netic  "dual,"  conducting  patches,  would  be  a  capacitive  screen.)  It  functions  as  a  high-pass 
filter,  reflecting  most  of  the  energy  below  the  frequency  at  which  the  aperture  spacing  is  one 
wavelength.  Combining  two  or  more  screens  in  series  sharpens  the  frequency  rejection,  by 
forming  a  cascaded,  two-pole  filter. 

Lee  et.  al.  solved  for  the  transmission  characteristics  of  these  screens  using  the  mode 
matching  method  [19].  This  method  expands  the  fields  in  a  series  of  plane  wave  modes  (Flo- 
Iquet  modes)  [20].  It  is  a  form  of  moment  method  in  the  sense  that  it  is  a  matrix  solution  to  a 
frequency-domain  integral  equation.  The  expansion  and  testing  functions  are  the  plane  wave 
modes  themselves  (entire  domain  functions). 

Figure  8  shows  the  representation  of  one  unit  cell  of  this  structure  by  a  tetrahedron  mesh. 
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The  shaded  areas  indicate  those  cell  faces  that  are  defined  as  perfect  conductor.  The  mesh 
density  is  approximately  13  cells  per  wavelength  at  the  highest  frequency,  which  corresponded 
to  a  unit  cell  width  and  height  of  one  wavelength.  For  this  test  case,  the  spacing  between 
screens  was  d=.2a  and  the  aperture  size  was  a=b=Jd^=Jdy. 

Figure  9  shows  the  calculated  power  transmission  through  the  double  screen.  The 
agreement  between  the  HFEM  code  and  the  moment  method  calculations  from  [19]  is  excellent. 
This  provided  an  initial  validation  of  several  features  of  the  theory,  such  as  the  periodic  radiation 
condition,  the  side-wall  periodicity  conditions,  and  the  interior  finite  element  calculations.  In 
short,  these  results  could  not  have  been  obtained  if  there  were  any  errors  in  the  theory  or 
implementation  of  the  matrix  equation  derived  in  Chapter  3.  It  also  verifies  the  correctness  of 
the  formulation  for  the  forcing  function  due  to  the  incident  field,  and  the  subsequent  calculation 


dx  /  A, 


Figure  9.  Comparison  of  HFEM  and  Method  of  Moments  [19]  Calculations  for 
Transmission  through  Two-Layer  Inductive  Screen:  d^=dy,  a—b^.ld^,  d=.2a 
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of  transmissivity  from  the  electric  field  solution. 

This  test  case  also  illustrates  that  the  FEM  does  not  have  any  difficulty  with  sharp 
conducting  edges,  and  does  not  require  extraordinarily  high  grid  resolution  near  such  edges.  It 
also  serves  to  illustrate  that  the  outer  surfaces  of  the  material  structure,  conducting  surfaces  in 
particular,  may  coincide  with  the  radiation  boundary.  Hence,  it  is  not  necessary  to  add  any  grid 
cells  in  free  space  above  or  below  the  structure. 

So2o  Wire  Grid  Array 

As  mentioned  earlier,  wire  grids  can  be  used  as  angular  filters  to  reduce  the  sidelobes 
of  an  antenna.  In  addition,  linearly  polarized  reflector  antennas  are  often  constructed  from 
parallel  wires  or  bars  in  order  to  reduce  weight  and  wind  resistance.  To  properly  design  such 
an  antenna,  it  is  necessary  to  know  how  closely  the  wires  must  be  spaced  to  provide  good 
reflection  at  all  incidence  angles. 

Decker  performed  measurements  of  transmission  through  a  wire  grid  as  a  function  of 
angle,  with  the  incident  electric  field  polarized  parallel  to  the  wires  [21],  He  compared  those 
measurements  to  analytical  formulas  from  Wait  [22]. 

To  use  the  hybrid  finite  element  method  for  this  calculation,  the  grid  shown  in  Figure 
10  was  used.  The  short  section  of  wire  is  represented  as  a  void  in  the  tetrahedron  mesh,  indicat¬ 
ed  by  the  shaded  area.  The  grid  is  fairly  wide  in  the  x  dimension,  corresponding  to'the  spacing 
between  wires  (1.5  cm).  The  wire  radius  is  1  mm. 

Figure  1 1  compares  the  HFEM  results  with  those  given  by  Decker,  both  measured  and 
calculated  (f=9.6  GHz).  The  HFEM  calculations  are  somewhat  closer  to  the  measured  data, 
most  notably  at  normal  incidence.  The  discrepancy  between  the  measurement  and  both  sets  of 
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Figure  10.  Tetrahedron  Mesh  for  Wire  Grid  Array 
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Figure  11.  Parallel-polarized  Transmission  through  Wire  Grid  Array  vs.  Scan  Angle 


calculations  at  wide  angles  is  most  likely  due  to  the  fact  that  the  test  article  was  fairly  small  in 
both  width  and  length,  allowing  for  diffractions  from  its  edges. 

It  is,  of  course,  quite  inefficient  to  use  a  three-dimensional  code  for  this  inherently  two- 
dimension  problem.  However,  it  serves  to  illustrate  that  the  HFEM  code  correctly  models 
situations  in  which  infinite  conductors  are  represented  as  objects  projecting  through  opposite  unit 
cell  walls.  A  later  section  presents  a  similar  case  in  which  wires  (or  bars)  continue  through 
opposite  unit  cell  faces  in  both  directions,  representing  a  grid  of  crossed  wires. 


Artificial  dielectrics  are  synthetic  delay  media  made  up  of  dielectric  or  metallic  obstacles 
imbedded  in  a  natural-dielectric  medium,  or  binder.  The  obstacles  are  usually  arranged  in  a 
regular  lattice  and  are  individually  much  smaller  than  the  wavelength  at  which  the  structure  is 
intended  to  function.  Artificial  dielectrics  have  the  advantage  that  they  may  have  refractive 
indices  that  are  much  higher,  for  a  given  weight,  than  natural  dielectrics. 

One  such  medium  is  an  array  of  metal  spheres.  When  arranged  in  a  cubic  lattice,  the 
effective  permittivity,  permeability,  and  index  of  refraction  may  be  estimated  analytically  using 
conformal  transformations.  Lewin  derived  the  following  expressions  [23]: 
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where  is  the  relative  permittivity  of  the  binder  material,  is  the  wavelength  in  the  binder, 
D  is  the  sphere  diameter,  and  V  is  the  unit  cell  volume.  These  formulas  are  expected  to  be  valid 
for  D/Xj  <  .2  [4],  From  (78)  and  (79),  the  effective  index  of  refraction  is 
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The  phase  delay  (in  deg.)  through  a  layer  of  thickness  t  is 
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Figure  12  shows  a  tetrahedron  mesh  of  the  air  space  surrounding  a  metal  sphere.  The 
sphere  diameter  is  2  m,  and  the  cube  is  4  m  on  each  side,  so  that  F  =  64  m^  .  Since  the  sphere 
is  an  impenetrable  conductor,  it  is  represented  as  a  void  in  the  mesh. 

The  formulas  (78)-(82)  apply  to  a  uniform,  infinite  medium.  However,  the  calculations 
are  for  a  structure  with  finite  thickness.  Therefore,  it  will  include  edge  effects  at  the  surface 


Figure  12.  Tetrahedron  Mesh  for  a  Conducting  Sphere  in  a  Cubic  Lattice 
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There  is  a  direct  analogy  between  the  quantum  states  in  a  crystal  and  microwave  frequencies  in 
a  (much  larger)  scaled  lattice  of  dielectric  obstacles  or  voids  in  a  solid  dielectric  [25].  The 
dielectric  structure  will  have  a  frequency  stop  band  corresponding  to  the  photonic  band  gap. 

An  example  bandgap  structure  is  an  array  of  layered  dielectric  rods.  Figure  14  shows 
the  cross  sectional  geometry  of  one  rod  in  an  array  whose  transmissivity  was  measured  by  Kelly 
et.  al.  [26].  The  array  was  up  to  ten  layers  of  these  rods  in  a  12.7  mm  x  12.7  mm  square 
lattice.  The  acrylic  tubes  have  relative  permittivity  ej.=2.55  and  inner  and  outer  radii  of  b  = 
3.175  mm  and  c  =  4.763  mm.  Each  may  have  a  pyrex  core  with  radius  a  =  3.0  mm  and 
e^=4.2.  There  is  a  small  air  gap  between  the  pyrex  core  and  the  acrylic  sheath.  In  Kelly’s 
measurements,  seven  of  the  rows  had  pyrex  cores  and  the  last  three  had  the  acrylic  tubes  only. 

For  purposes  of  the  finite  element  calculations,  the  air  gap  is  ignored  and  the  acrylic  tube 
is  assumed  to  have  an  inner  diameter  of  6.0  mm.  The  tetrahedron  mesh  for  a  single  rod  is 
shown  exploded  in  Figure  15a.  Note  that  since  pyrex  has  a  refractive  index  of  approximately 
2,  the  mesh  in  that  material  is  twice  as  dense  as  in  the  air  region.  Figure  15b  shows  the  unit 
cell  mesh  for  a  four-layer  structure.  Only  the  very  narrow  end  of  the  mesh  is  subject  to  the 
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Figure  14.  Layered  Dielectric  Rod  used  to  Simulate  Photonic  Band  Structure 
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Figure  16.  Comparison  of  Hybrid  Finite  Element  Method  Calculations 
and  Measurements  for  10-layer  Dielectric  Rod  Array 

the  most  interesting  photonic  band  structures  are  crystal  lattices,  which  cannot  be  represented 
as  2-dimensional  structures.  The  hybrid  finite  element  method  clearly  has  the  capability  of 
dealing  with  general  3-dimensional  bandgap  structures,  and  the  potential  to  supplant  tedious  and 
expensive  measurements. 

5.5.  Reinforced  Concrete 

Reinforced  concrete  walls  are  usually  made  up  of  regular  lattices  of  steel  reinforcing  bars 
in  the  center  of  the  concrete  layer.  Their  level  of  protection  from  radio  frequency  energy,  or 
shielding  effectiveness,  has  much  more  to  with  the  reinforcing  bars  than  the  concrete  itself. 
There  is  tremendous  variation  in  the  electrical  properties  of  concrete,  depending  mainly  on  its 
water  content,  but  it  does  not  usually  attenuate  RF  energy  very  well  by  itself.  Hence,  predicting 
the  shielding  effectiveness  of  reinforced  structures  depends  on  being  able  to  predict  the  trans- 
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Evidently,  when  the  wave  is  polarized  parallel  to  either  set  of  bars,  the  transmissivity  is  the 
same  whether  or  not  they  are  electrically  connected.  However,  for  oblique  incidence,  the 
situation  is  quite  different,  as  Figure  19  attests.  Here,  the  wave  is  incident  from 

transverse  electric  polarization  (electric  field  polarized  45°  to  both  sets 
of  bars).  When  the  bars  are  connected,  circulating  currents  are  induced  that  create  a  resonance 
condition  near  1.75  GHz,  with  zero  co-polarized  transmission.  (There  is  a  small  amount  of 
power  transmitted  in  the  orthogonal  polarization.)  This  transmission  null  will  shift  in  frequency 
as  the  incidence  angle  changes.  Both  curves  indicate  a  drop  in  transmissivity  above  2.25  GHz, 
but  that  is  due  to  the  fact  that  the  transmitted  power  is  split  between  the  main  lobe  and  a^grating 
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Figure  17,  Cutaway  of  One  Half  of  Tetrahedron  Mesh  for  Reinforced  Concrete  Unit  Cell 

lobe.  The  grating  lobe  power  is  not  accounted  for  in  this  calculation.  The  total  transmissivity, 
including  the  grating  lobe,  will  be  near  unity  above  2.25  GHz. 

5.6.  Crossed  Dipole  Frequency  Selective  Surface 

The  crossed  dipole  FSS  illustrated  in  Figure  2  is  the  only  test  case  attempted  so  far  for 
which  satisfactory  results  have  not  yet  been  obtained.  The  crossed  dipole  array  was  an  early 
attempt  at  a  "dichroic"  surface  that  would  reflect  a  single  frequency  and  pass  all  others.  The 
intended  use  was  a  subreflector  for  a  large  parabolic  dish.  The  subreflector  would  pass  energy 
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from  a  feed  located  at  the  dish’s  focus,  while  reflecting  that  from  a  feed  located  in  the  center 
of  the  dish.  Unfortunately,  the  crossed  dipole  element  was  deficient  because  at  off  broadside 
scan  angles  it  has  an  anti-resonance,  sometimes  called  a  "Wood’s  anomaly"  [28].  Figure  20 
shows  the  element  and  lattice  dimensions  and  the  unit  cell  outline.  Figure  21  compares  the 
HFEM  calculations  with  method  of  moment  results  from  Cwik  &  Mittra  [29].  The  effect 
observed  when  attempting  this  problem  with  HFEM  is  a  shift  in  the  resonant  frequency  and  a 
broadening  of  the  frequency  response.  The  grid  used  for  the  HFEM  results  is  shown  in  Figure 
22,  which  includes  two  air  layers  on  each  side  of  the  dipole  elements.  The  grid  had  a  mesh 
edge  length  less  than  X/30  in  the  plane  of  the  dipoles  at  the  frequency  where  the  anti-resonance 
occurs.  Even  though  that  would  ordinarily  be  an  excessively  fine  mesh,  it  is  evidently  not  fine 
enough  in  this  case.  Accurate  results  for  such  cases  may  require  the  development  of  special 
singularity  finite  elements.  Although  such  elements  are  common  in  mechanical  engineering 
applications,  they  are  only  beginning  to  see  development  for  electromagnetic  problems. 

5.7.  Storage,  Execution  Time  and  Accuracy  Issues 

Table  I  shows  the  sizes  of  each  of  the  test  cases  in  terms  of  the  mesh  sizes  and  storage 
required.  In  all  cases,  the  storage  needed  for  the  matrix,  both  sparse  and  dense  parts,  is  less 
than  4  Mb  (Megabytes)  using  single-precision  complex  representation  (8  bytes  per  number,  or 
"word").  The  total  number  of  mesh  edges  is  almost  always  less  than  1.5  times  the  number  of 
mesh  cells.  For  the  sparse  matrix,  there  are  typically  15  or  less  entries  per  matrix  row.  The 
number  of  edges  in  the  two  radiation  boundaries  is  always  less  than  where  f  is  the 

mesh  edge  length.  An  exception  is  when  part  of  the  radiation  boundary  is  filled  with  conducting 
surfaces,  as  in  the  case  of  the  inductive  grid  (see  Fig.  8).  In  summary,  the  storage  required  for 
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Figure  22.  Tetrahedron  Mesh  for  Crossed  Dipole  Frequency  Selective  Surface 
Test  Case;  Two  Air  Layers  Surround  the  Conducting  Elements 

any  problem  may  be  estimated  a  priori  as 

Matrix  Storage  {words)  <  (15)(1.5)  (83) 

where  is  the  number  of  mesh  cells. 

Table  II  summarizes  the  execution  time  for  each  case  using  a  SparcStation™-20  comput¬ 
er.  These  are  averages  over  all  of  the  frequencies  or  angles  for  which  results  were  shown  in 
the  previous  sections.  Each  new  angle  or  frequency  requires  a  separate  matrix  fill  and  solve. 
It  is  immediately  apparent  that  the  biconjugate  gradient  solver  has  a  distinct  advantage,  converg- 
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ing  5-12  times  faster  than  the  conjugate  gradient  solver.  Unfortunately,  the  BCGM  did  not 
converge  at  all  for  several  frequencies  of  the  bandgap  structure  test  case.  It  also  failed  to 
converge  for  any  frequency  of  the  crossed  dipole  FSS  test  case. 

Figure  23  shows  the  convergence  vs.  frequency  for  the  bandgap  structure.  CGM’s 
variation  is  fairly  regular,  suggesting  a  systematic  dependence  on  some  feature  of  the  problem 
that  has  not  yet  been  identified.  BCGM’s  variation  with  frequency  is  more  random.  The  circles 
along  the  top  of  the  graph  mark  those  frequencies  where  BCGM  failed  to  converge. 

These  execution  statistics  show  that  the  HFEM  technique  can  solve  many  practical 
problems  in  a  reasonable  time  on  a  typical  engineering  workstation.  A  final  important  note  is 
that  its  efficiency  is  much  poorer  than  moment  methods.  Its  use  is  justified  when  the  problem 
geometry  includes  features  that  MoM  cannot  model. 


Figure  23.  Number  of  Iterations  vs.  Frequency  for  Bandgap  Structure  Test  Case 
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6o  CONCLUSIONS  AND  RECOMMENDATIONS 

The  objective  of  this  project  in  computational  electromagnetics  was  to  develop  a  method 
for  predicting  reflection  and  transmission  from  general  periodic  structures.  The  hybrid  finite 
element  method  was  selected  because  of  its  inherent  ability  to  deal  with  inhomogeneous  dielec¬ 
trics  and  arbitrarily-oriented  conductors.  This  report  has  discussed  the  theoretical  basis  in  terms 
of  the  electric  field  wave  equation  with  appropriate  periodic  boundary  conditions,  its 
into  a  numerical  problem  and  matrix  solution,  its  implementation  in  a  general-purpose  computer 
program,  and  finally,  its  results  for  a  wide  variety  of  validation  cases.  Those  results  show  that 
the  method  is  effective  and  reliable,  although  its  efficiency  relative  to  other  methods  is  problem- 
dependent.  The  variety  of  validation  cases  accomplished  with  a  single  computer  code  demon¬ 
strates  that  HFEM  is  the  most  versatile  tool  to  date  for  periodic  structure  analysis. 

Recommendations  for  further  work  fall  into  two  categories:  exploiting  the  new  code  for 
problems  of  immediate  interest  to  high  power  microwave  research;  and  improvements  to  extend 
the  code’s  range  of  usefulness  and  validity.  Besides  its  demonstrated  utility  for  reinforced 
concrete  walls,  the  code  may  also  be  used  to  predict  the  shielding  effectiveness  of  brick  walls, 
including  air  pockets,  perforated  metal  grids  such  as  those  used  for  equipment  enclosures,  and 
radomes  and  canopies  that  incorporate  metallic  grids  for  interference  suppression.  Artificial 
dielectrics  may  also  be  interesting  for  HPM  antenna  design,  possibly  replacing  heavier  natural 
dielectric  lenses.  But  since  the  artificial  dielectrics  cause  the  fields  passing  through  to  become 
concentrated  in  between  conducting  elements,  they  may  experience  dielectric  breakdown  at  much 
lower  field  levels.  HFEM  provides  a  unique  tool  for  examing  the  electromagnetic  field  distribu¬ 
tion  in  candidate  structures  and  determining  which  are  likely  to  be  most  useful, 

RF  shielding  by  composite  materials  made  of  graphite  fibers  embedded  in  a  dielectric 
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material  such  as  epoxy  will  become  interesting  in  the  future.  For  modeling  those  structures,  it 
is  more  efficient  to  represent  the  graphite  fibers  as  conducting  filaments  with  a  surface  resistivi¬ 
ty.  Hence,  a  proposed  extension  to  the  computer  code  is  to  add  conductors  with  finite  conduc¬ 
tivity.  Second,  it  was  observed  that  some  problems  experienced  very  slow  convergence  using 
the  iterative  matrix  solvers.  Those  would  benefit  from  implementation  of  a  direct  matrix  solver 
using  LU  decomposition.  However,  that  solver  must  take  advantage  of  the  matrix  sparsity  by 
starting  with  a  row/column  reordering  such  as  the  "minimum  degree"  algorithm  [30]  to  minimize 
matrix  fill  in  as  the  decomposition  proceeds.  Third,  this  periodic  structure  code  is  one  of  a 
family  of  three  HFEM  codes,  the  other  two  of  which  perform  analyses  of  waveguide  devices 
[31]  and  phased  array  antennas  [6].  It  is  both  possible  and  desirable  to  combine  those  capabili¬ 
ties  in  a  single  code.  Finally,  the  crossed-dipole  FSS  test  case  illustrated  a  shortcoming  of  the 
present  implementation  of  vector  finite  element  solutions  in  regions  that  include  field  singulari¬ 
ties.  As  an  alternative  to  very  fine  meshing,  which  is  computationally  expensive,  expansion 
functions  ("singularity  elements")  that  can  model  the  singular  behavior  are  preferred  and  should 
be  devloped. 
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